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Abstract 

This  paper  addresses  the  development  of  energy-based  models  and  model-based  control  designs  necessary 
to  achieve  present  and  projected  applications  involving  atomic  force  microscopy.  The  models  are  based  on  a 
combination  of  energy  analysis  at  the  mesoscopic  level  with  stochastic  homogenization  techniques  to  construct 
low-order  macroscopic  models.  Approximate  model  inverses  are  then  employed  as  filters  to  linearize  transducer 
responses  for  linear  robust  control  design. 
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1.  Introduction 

Developed  in  1986,  the  atomic  force  microscope  (AFM)  relies  on  interatomic  forces  between  a  cantilever 
tip  and  the  sample  to  obtain  ultra-high  resolution  surface  images  [4],  The  relatively  low  cost  of  the  devices 
and  the  fact  that  they  require  minimal  sample  preparation  has  made  the  AFM  a  standard  diagnostic  tool  in 
research  laboratories.  However,  several  present  and  projected  applications  make  requirements  on  the  technology 
that  current  AFM  designs  are  not  able  to  consistently  achieve.  These  limitations  are  primarily  due  to  the 
high  sample  rates  required  for  applications  such  as  real-time  product  diagnostics  or  monitoring  of  biological 
processes,  nanoelectromechanical  (NEMS)  applications,  and  employment  of  AFM  technologies  for  spintronics. 
Real-time  product  diagnostics  include  analysis  of  contact  lenses  to  detect  defects  or  protein  deposits  and  screening 
of  semiconductor  chips  to  maintain  quality  control.  The  real-time  monitoring  of  biological  processes  has  the 
potential  for  leading  to  treatment  policies  for  ailments  such  as  osteoporosis  as  well  as  the  potential  for  quantifying 
fundamental  biological  phenomena  such  as  protein  unfolding.  Within  the  NEMS  regime,  the  repulsive  forces 
utilized  in  atomic  force  microscopy  also  lead  to  the  potential  for  nanoconstruction  using  the  cantilever  as  an 
actuator.  Finally,  the  extreme  accuracy  provided  by  the  AFM  is  presently  being  combined  with  nuclear  magnetic 
resonance  microscope  (NMRM)  technology  to  investigate  the  detection  of  single  electron  spins  [3,  8,  16]. 

A  crucial  component  in  the  AFM  design  for  all  of  these  applications  is  the  piezoceramic  (PZT)-based  stage 
used  to  position  the  sample  as  depicted  in  Figure  1(a)  [2,  4,  11,  13,  14].  Nested  minor  loops  in  data  collected 
at  0.1  Hz  are  shown  in  Figure  1(b).  Whereas  PZT  actuators  provide  the  broadband  and  extremely  high  set 
point  capabilities  required  by  the  AFM  stages,  they  also  exhibit  frequency-dependent  hysteresis  as  illustrated 
in  Figure  2.  At  low  frequencies,  the  hysteresis  inherent  to  the  these  materials  can  be  accommodated  through 
PID  or  robust  control  designs  [1,  9].  However,  at  the  higher  frequencies  required  by  the  previously  outlined 
applications,  increasing  noise-to-data  ratios  and  diminishing  high-pass  characteristics  of  control  filters  preclude 
a  sole  reliance  on  feedback  laws  to  eliminate  hysteresis.  This  motivates  the  development  of  control  designs  that 
incorporate  and  approximately  compensate  for  hysteresis  through  model  inverses  employed  either  in  feedback  or 
feedforward  loops. 
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(a)  (b) 

Figure  1:  (a)  PZT-based  AFM  stage,  (b)  Nested  minor  loops  in  data  collected  at  0.1  Hz  on  an  AFM. 
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Figure  2:  Frequency  dependent  field-displacement  behavior  of  an  AFM.  Sample  rates  of  (a)  0.279  Hz,  (b)  1.12  Hz, 
(c)  5.58  Hz,  and  (d)  27.9  Hz. 
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Figure  3:  Gibbs  energy  G  and  hysteron  P  resulting  from  the  necessary  condition  fp  =  0  for  negligible  thermal 
activation. 


2.  Constitutive  Relations 

To  model  the  constitutive  behavior  of  the  piezoceramic  stacked  actuator,  the  stress-strain  relation  is  assumed 
to  be  linear.  However,  the  relation  between  the  applied  voltage  V  (or  the  applied  field  E)  and  the  polarization  P 
exhibits  nonlinearities  and  hysteresis.  The  actuator  is  also  assumed  to  be  biased  through  poling  so  that  the 
relation  between  P  and  the  strain  e  is  linear  for  the  considered  operating  conditions.  To  characterize  the 
hysteretic  E-P  behavior  at  the  domain  level,  a  Helmholtz  energy  relation  was  derived  in  [15]  using  statistical 
mechanics  principles  under  the  assumption  that  dipoles  are  either  aligned  with  the  field  or  diametrically  opposed 
to  it.  This  model  is  appropriate  for  a  single  crystal  with  uniform  effective  fields.  To  construct  a  macroscopic 
model  for  a  polycrystalline  material  with  variable  effective  fields,  the  coercive  and  effective  field  values  are  then 
assumed  to  be  distributed  as  detailed  in  the  latter  half  of  this  section. 

Under  fixed  temperature  conditions  with  no  applied  stress  er,  it  is  illustrated  in  [15]  that  a  first  order  approx¬ 
imation  to  the  statistical  mechanics-based  Helmholtz  energy  is  the  piecewise  quadratic  relation 

(±r,(P  +  PR)2  P<~Pi 

Hp)  =  l  HP  -  PR )2  P  >  Pi  (1) 

{Hpi  -  Pr)  ~  Pr)  \p\<Pr 

As  shown  in  Figure  3,  Pi  is  the  positive  inflection  point  and  PR  is  the  value  of  P  at  which  the  positive  local 
minimum  of  ip  occurs.  The  parameter  77  is  the  reciprocal  of  the  slope  of  the  E-P  relation  after  switching  occurs. 
This  fact  can  be  used  to  establish  an  initial  parameter  value  for  77  when  modeling  a  specific  data  set. 

In  the  case  of  no  applied  stresses  er,  the  Gibbs  energy  can  be  formulated  as 

G{E,  P)  =  ip  —  EP  (2) 

where  the  second  term  represents  the  electrostatic  energy  due  the  applied  field  E.  In  order  to  include  ferroelastic 
coupling,  we  utilize  the  extended  Helmholtz  relation 

H(P,  e)  =  HP)  +  \yPz2  ~  YPieP.  (3) 


3 


The  Gibbs  energy  is  then  given  by 


G(E,  P,  e)  =  ijj(P)  +  ^ Ype 2  -  Yp^eP  -  EP  -  oe 


(4) 


where  ae  incorporates  the  elastic  energy.  Note  that  Yp  is  the  Young’s  modulus  for  a  constant  polarization  and 
7  is  a  ferroelastic  coupling  coefficient. 

In  the  case  of  negligible  thermal  activation,  the  local  average  polarization  P  is  determined  from  the  necessary 
conditions 


d  G 
dP 


=  0 


d2G 

dp2 


>0. 


Applying  these  conditions  to  (2)  yields  a  piecewise  linear  E-P  characterization 


(5) 


r[P(P;Pc,0](0),  r(t)  =  0 

[P(E;  Ec,£)](t)  =  <  |  -  Pr,  r(t)  ^  0  and  £(maxr(t))  =  -Ec  (6) 

[  f  +  Pr,  r(t)  ^  0  and  £(maxr(t))  =  Ec 

where 

ff  -Pr,  E(0)<-Ec 

[■ P(E ;  Pc,  0]  (0)  =  <  £,  -Pc  <  p(o)  <  Pc  (7) 

If  +  Pr,  P(0)>Pc 

and  the  transition  times  are 

r{t)  =  {t£  (0,  Tf]  I  Et  =  -Pc  or  Et  =  Ec}.  (8) 

However,  if  thermal  activation  is  significant,  the  dipoles  can  achieve  the  thermal  energy  required  to  switch 
in  advance  of  the  minimum  Gibbs  energy  so  the  relative  thermal  and  Gibbs  energy  must  be  balanced  through 
Boltzmann  principles.  The  probability  density  for  achieving  an  energy  level  G  is  then  given  by 


H(G)  =  Ce~GV/kT 


(9) 


where  k  is  Boltzmann’s  constant,  V  is  a  reference  volume  and  C  is  a  constant  that  is  selected  so  that  when  /i(G) 
is  integrated  over  all  possible  dipole  orientations,  a  probability  of  1  is  achieved.  If  we  let  2<r  be  the  separation 
between  possible  polarization  states  around  P0,  the  probabilities  of  reaching  a  polarization  state  having  sufficient 
energy  to  switch  orientations  are  given  by 


rPo+°  G{E,P)V/kTdp 

_  JPp-<y _ 

e-G(E,P)V/kTdp 

J  fQ  —  <7 


cPo+a  —G(E,P)V/kT  dp 
r_  _  JPo-v _ 

+  fPo+<T  e-G(E,P)V/kTdp  ' 
J  —  oo 


(10) 


The  likelihoods  of  reaching  the  required  energy  and  thus  of  the  dipoles  switching  from  a  positive  to  a  negative 
orientation  and  conversely  are  then 


1  1 

P+-  =  ~r+-  ,  P-+  =  ~r-+ 

T  T 


(11) 


where  r  is  the  relaxation  time.  The  fractions  of  dipoles  in  each  orientation  evolve  according  to  the  ordinary 
differential  equations 


dx+  dx _ 

——  =  -p+_x++p_+x-  ,  — —  =  - p-+x _  +p+_x+. 


The  expected  polarizations  due  to  positively  and  negatively  oriented  dipoles  are 

coo  rPo— <r 

(P+)  =  /  Pp{G)dP  ,  (P+)  =  /  Pp(G)dP 

J  Pq-\-<j  J —  oo 


(12) 


(13) 
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so  the  evaluation  of  C  yields 

f“  pe-G(E,P,T)V/k.Tdp  rP b-°  pe-G(E,P,T)V/kTdp 

+  /  Jpo+ae~G(E^T)V/kTdP  ’  V  e-G(E,P,T)V/kTdp  ' 

For  a  single  crystal  with  uniform  effective  field,  the  local  average  polarization  is  subsequently 

P  =  x+{P+)+x-(P _>. 

In  the  manner  detailed  in  [15],  the  evaluation  of  the  integrals  in  (10)  and  (14)  can  be  simplified  through  approx¬ 
imations  employing  the  inflection  points  ±Pj  rather  than  the  unstable  equilibrium  Pq. 

Both  of  the  relations  (6)  and  (15)  are  valid  only  for  homogeneous  single  crystal  materials  with  uniform  effective 
fields.  To  account  for  nonuniformity  and  inhomogeneities  in  the  materials,  local  coercive  and  effective  fields  are 
assumed  to  be  manifestations  of  underlying  distributions  rather  than  constants.  The  macroscopic  polarization 
model  is  then  given  by 

/>00  POO 

P(E)  =  /  /  P(E  +  Ee,Ec,^u1(Ec)p2(Ee)dEedEc  (16) 

J  0  J  —  oo 

where  zq  and  zq  are  appropriate  densities.  Motivated  by  choices  in  the  magnetics  literature,  zq  and  zq  were 
respectively  designated  to  be  lognormal  and  normal  densities  in  [15].  However,  the  fact  that  neither  of  these 
choices  is  based  on  energy  considerations,  motivates  the  consideration  of  general  densities. 

Based  on  physical  arguments,  we  may  assume  that  the  general  densities  decay  to  zero.  Therefore,  the  double 
integral  can  be  evaluated  by  truncating  the  domains  and  using  composite  Gaussian  quadrature.  In  this  case, 
discretization  of  (16)  gives 


(14) 

(15) 


Nt  Nj 

P{E)  =  y  P(E  +  Eei ,  ECi ,  &) V\ (Ee. ) u2 (ECi)vjWj 

i= 1  3= 1 

(17) 

where  Ee.,Ec.  are  the  abscissas  and  Vi,Wj  are  the  weights.  Here,  v\(Ee.)  and  V2{EC .) 
be  determined.  Details  regarding  the  identification  of  zq  and  zq  can  be  found  in  [10]. 

are  parameter  values  to 

Finally,  the  equilibrium  condition 

^=0 

OS 

(18) 

yields  the  elastic  constitutive  relation 

a  =  Ype  -  Yp^P. 

(19) 

This  relation  along  with  the  nonlinear  polarization  relation  (16)  quantifies  the  constitutive  behavior  for  the 
piezoceramic  materials  employed  in  the  AFM. 


3.  Model  Inverse 

To  construct  a  model  inverse  to  be  employed  as  filter  in  the  robust  control  design,  consider  first  an  imple¬ 
mentation  of  the  forward  model  (17)  described  in  [15]  and  summarized  below.  Formulate  the  local  polarization 
(6)  as 

P=-+PrA  (20) 

V 

where  A  is  an  V,;  x  Nj  matrix  whose  elements  are  either  1  or  —1.  The  ij th  entry  in  A  indicates  whether  the  jth 
effective  field  value  Eej  has  crossed  the  ith  coercive  field  value  ECi  and  thus  whether  the  associated  polarization 
value  is  on  the  upper  or  lower  branch  of  the  hysteron.  Now  rewrite  A  as  a  vector  Avec  with  the  order  of  the 
entries  determined  by  order  in  which  the  entries  in  A  change  as  the  field  E  is  increased  or  decreased.  Thus  Avec 
can  be  updated  one  entry  at  a  time  as  the  field  changes. 


5 


A  model  for  quantifying  the  inverse  map  between  P  and  E  is  now  outlined  in  Algorithm  1. 


Algorithm  1. 


Specify  Eprev 


P  P 

1  prev ?  1  n 


temp  —  *  prev 
while  Ptemp  Pprev 

Etemp  =  value  of  E  at  which  the  next  entry  in 
Update  next  entry  in  Avec 
Compute  change  in  P  and  update  P 
Update  Ptemp 
end 

Compute  Enew 


A,; 


changes 


4.  Stacked  Actuator  Model 


In  addition  to  hysteresis,  the  dynamics  of  the  actuator  must  be  incorporated.  We  assume  that  the  rod  has 
cross-sectional  area  A ,  length  £,  density  p  and  Young’s  modulus  Yp .  Let  cp  be  the  Kelvin- Voight  damping 
parameter  and  7  be  the  piezoelectric  coupling  coefficient.  In  the  present  stage  design,  depicted  in  Figure  1(a), 
one  end  of  the  actuator  is  fixed,  while  the  attachment  at  the  other  end  can  be  modeled  as  a  damped  spring-mass 
system.  For  this  end,  let  M L  be  the  mass,  k l  be  the  stiffness,  and  cl  be  the  damping  coefficient.  Force  balancing 
along  the  actuator  then  yields  the  relation 


d2u  dJ\f 


(21) 


where  the  resultant  M  =  fA  crdA  is  given  by 


(22) 

Note  that  u  denotes  displacement  in  the  longitudinal  ^-direction  and  that  the  relation  strain  s  =  ^  is  used  in 
obtaining  (22).  The  boundary  conditions  are  u(t,  0)  =  0  at  the  fixed  end  and 


A f(t,£) 


,  ,  „  du ,  , ,  d2u , 

-kLu(t,£)  -  cL  —  (t,£)  - 


(23) 


at  the  moving  end.  Initial  conditions  are  given  by  rt(0,a;)  =  =  0.  The  polarization  P(E)  is  specified  by  (16). 

An  approximate  solution  to  (21)  is  found  by  first  deriving  a  weak  form 


ipdx  + 


YpA—  +  cpA— 
dx  dxdt 


^-dx  =  -  [  Yp A(3P  {E)  -7r-dx 
c)x  Jq  c)x 


(24) 


and  then  discretizing  in  space  with  linear  finite  elements  and  in  time  with  finite  difference  techniques.  Details  of 
the  derivation  of  the  weak  form  and  the  construction  of  the  finite  element  and  finite  difference  equations  can  be 
found  in  [13]. 

In  order  to  determine  the  displacement  of  the  rod  at  the  end  x  =  £,  the  rod  model  can  be  approximated  by 
a  lumped  spring-mass  model  since  the  cross-section  is  small  compared  to  the  length  and  the  electric  field  E  is 
approximately  uniform  throughout  the  actuator.  Numerical  simulations  and  physical  experiments  demonstrate 
that  the  resulting  ODE  model  is  a  good  approximation  to  the  PDE  model  for  specifying  tip  displacement.  The 
ODE  model  has  the  advantage  of  faster  computational  times  compared  to  the  discretized  PDE. 
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Figure  4:  Characterization  of  AFM  field-displacement  behavior  at  0.1  Hz. 


5.  Model  Validation 

To  demonstrate  the  accuracy  and  efficiency  of  the  hysteresis  model,  we  consider  the  characterization  of 
nested  minor  loops  collected  at  0.1  Hz  as  well  as  E-P  behavior  at  frequencies  ranging  from  0.28  Hz  to  27.9  Hz. 
In  both  cases,  displacements  were  computed  using  the  ODE  model  constructed  using  the  stress  relation  (19)  and 
hysteresis  model  (16)  with  general  densities  v\  and  identified  using  the  techniques  detailed  in  [10,  12]. 

Figure  4  illustrates  the  capability  of  the  model  to  characterize  nested,  biased  minor  loop  behavior  and  Figure  5 
demonstrates  the  characterization  of  frequency-dependent  dynamics.  The  latter  involves  the  quantification  of 
both  thermal  relaxation  and  inertial  effects  as  illustrated  by  the  change  in  sign  of  the  slope  following  field 
reversal.  Further  details  demonstrating  properties  of  the  model  for  characterizing  hysteresis  in  various  PZT 
compounds  can  be  found  in  [10,  11,  12,  15]. 


6.  Robust  Control  Design 

To  construct  a  robust  control  design  that  accommodates  the  hysteresis  and  constitutive  nonlinearities,  the 
model  can  be  inverted  in  the  manner  outlined  in  Section  3  and  employed  as  a  filter.  Discretization  errors  combined 
with  modeling  errors  will  prevent  complete  elimination  of  error  by  the  filter,  but  the  disturbance  d  associated 
with  the  actuator  can  be  significantly  reduced  in  this  manner.  We  will  construct  TL2  and  7foo  algorithms  that 
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Figure  5:  Characterization  of  AFM  field-displacement  behavior  with  sample  rates  of  (a)  0.279  Hz,  (b)  1.12  Hz, 
(c)  5.58  Hz,  and  (d)  27.9  Hz. 


utilize  the  model  inverse.  After  a  discussion  of  the  system  representation,  numerical  examples  are  given  to 
demonstrate  the  performance  of  the  control  designs. 

The  control  design  is  analogous  to  that  described  in  [5,  6,  7].  The  physical  control  system  is  the  AFM  and  the 
control  objective  is  to  track  a  reference  trajectory  r  which  represents  the  displacement  of  the  actuator.  The  system 
representation  in  shown  in  Figure  6.  The  plant  P  represents  the  ODE  model  for  the  AFM  rod.  The  sensor  noise 
in  the  measurement  of  y  is  separated  into  two  components,  s,  taken  to  be  60  Hz  noise  to  simulate  electromagnetic 
disturbances,  and  n ,  taken  to  be  high  frequency  noise  which  can  be  attributed  to  the  sensing  device  or  external 
disturbances.  The  disturbance  d  represents  errors  in  the  plant  input  caused  by  either  unattenuated  hysteresis 
and  constitutive  nonlinearities  or  by  discretization  errors  from  the  inverse  filter.  The  output  signal  e  represents 
the  weighted  tracking  error  and  u  is  the  weighted  output  of  the  controller  K .  The  weighting  functions  Wr,  We, 
Wv ,  Wd,  Ws,  and  Wn  are  selected  to  achieve  the  best  performance  from  the  controller.  They  are  chosen  based 
on  known  information  about  the  corresponding  signals. 

To  describe  the  open  loop  system,  maps  from  the  inputs  r,  d ,  s,  and  n  to  the  outputs  e,  e  and  u  are  designated 
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Figure  6: 


System  representation  including  input  disturbance  d  and  sensor  noise  s  and  n  in  the  transducer. 


as 


e  =  Wr[r\-(P[Wd[d\+u}+Wn[n}  +  Ws[s\ 
=  Wr [r]  -  P[Wd [d\]  -Wn-Ws [s]  -  P[u] 
e  =  We[e] 
u  =  Wu[u]. 


The  transfer  function  matrix  G  from  the  inputs  r,  d ,  n,  s  and  u  to  the  outputs  e,  e  and  u  is  then  specified  by 


WeWr  -WePWd  -wewn  -wews 
0  0  0  0 
Wr  -PWd  -Wn  -Ws 


-WeP 

wu 

p 


(25) 


The  linear  fractional  transformation  system  representation  is  shown  in  Figure  7  with  details  provided  in  [7,  17]. 

The  choice  of  weighting  functions  is  an  important  part  of  the  control  design.  To  filter  the  60  Hz  sensor  noise  s, 
a  sixth-order  passband  Chebyshev  filter  with  a  range  of  55  Hz  to  65  Hz  is  used  for  Ws.  The  frequency  response 
of  Ws  is  shown  in  Figure  8(a).  The  choice  for  Wn  is  a  second-order  highpass  Butterworth  filter  with  a  cutoff  of 
150  Hz  which  accommodates  the  high  frequency  sensor  noise  n.  The  weight  Wr  is  based  on  which  components 
of  the  reference  signal  are  most  important  during  tracking.  For  example,  considering  the  0.25  Hz  sinusoidal  part 
of  the  reference  signal  to  be  most  important,  a  sixth-order  passband  Chebyshev  filter  with  a  range  of  0.125  Hz 
to  0.375  Hz  was  employed.  To  determine  the  weighting  filter  Wd,  the  reference  signal  was  scaled  to  have  the 
same  order  of  magnitude  as  polarization.  This  signal  was  then  fed  into  either  the  inverse  filter  or  a  linear  filter 
to  determine  the  scale  from  polarization  and  field.  The  result  was  then  fed  to  the  forward  model.  This  process 
is  depicted  in  Figure  9  for  the  two  choices  of  disturbance.  In  both  cases,  the  power  spectrum  of  the  output  d 
had  only  low  frequency  components.  Thus  a  second-order  lowpass  Butterworth  filter  with  a  cutoff  of  10  Hz  was 
used  for  Wd.  The  frequency  response  of  Wd  in  the  case  when  the  linear  filter  was  used  is  shown  in  Figure  8(b). 
The  weighting  function  on  the  error  signal  e  was  chosen  to  be  We  =  with  qe  =  1  x  101  and  ee  =  1  x  1CA8. 
An  integrator  was  chosen  so  that  the  error  would  not  achieve  steady  state  at  a  nonzero  value.  Also,  the  pole 
was  shifted  slightly  off  zero  to  ensure  that  the  controller  was  realizable.  Finally,  the  weighting  function  on  the 
controller  output  was  chosen  to  be  Wu  =  5  x  10-6. 


r 

d 

s 

n 


=  w 


Figure  7:  Linear  fractional  transformation  representation  (LFT)  of  the  transducer  model. 
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(a)  (b) 

Figure  8:  (a)  Frequency  response  of  the  passband  filter  Ws  (b)  Frequency  response  of  W,i  for  the  disturbance  d 
due  to  scaled  but  uncompensated  hysteresis  and  nonlinearities. 

The  transfer  function  matrix  (25)  gives  a  representation  for  the  transducer  system  with  weighting  filters 
included.  In  order  to  formulate  the  control  laws  for  computing  the  gains  K,  7Y2  and  Hoo  norms  of  the  closed 
loop  system  representation  T  are  minimized.  These  are 

\\T\\l  =  ip  /  trace  [T*(ju)T(ju)\  dui  ,  HTH^  =  sup  d  [T(ju)]  (26) 

Zn  J-oo  U6R 

where  a[T(ju>)]  represent  the  maximal  singular  values  of  the  closed  loop  map  T .  Details  regarding  this  robust 
control  formulation  can  be  found  in  [7,  17]. 


7.  Numerical  Examples 

We  consider  the  two  possibilities  for  the  disturbance  depicted  in  Figure  9.  The  first,  which  is  due  to  scaled 
but  uncompensated  hysteresis,  yields  the  tracking  results  and  errors  shown  in  Figure  10.  In  the  second  case,  the 
disturbance  is  due  to  errors  existing  in  the  inverse  compensation  procedure  used  to  approximately  linearize  the 
transducer  response.  The  tracking  capabilities  and  errors  for  this  case  are  shown  in  Figure  11. 

A  comparison  between  Figures  10  and  11  illustrates  that  highly  accurate  tracking  is  obtained  in  both  cases, 
with  errors  less  than  2  (im  maintained  after  the  commencement  of  the  periodic  cycle.  The  equivalence  in  accuracy 
for  the  uncompensated  and  compensated  designs  is  attributed  to  the  low  level  of  hysteresis  present  in  this  low 
frequency  0.25  Hz  drive  regime  as  demonstrated  by  the  data  in  Figure  2.  The  present  investigation  focuses  on 
the  extension  of  these  control  designs  to  higher  frequencies  where  increasing  hysteresis  levels  necessitate  inverse 
compensation  as  demonstrated  for  analogous  magnetic  model  based  controllers  in  [7]. 

Very  similar  results  are  obtained  using  the  'Hco  design,  indicating  that  for  this  application,  the  construction 
of  models  and  control  filters  may  play  a  more  important  role  than  the  choice  of  robust  control  laws  for  high 
accuracy  tracking. 


Figure  9:  (a)  Disturbance  d  due  to  scaled  but  unmodeled  hysteresis  and  constitutive  nonlinearities,  (b)  Distur¬ 
bance  d  due  to  inverse  filtering  errors. 
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(a) 


Time 

(b) 


Figure  10:  TI2  design  with  sensor  noise  s  and  the  disturbance  d  due  to  inversion  error,  (a)  Reference  and 
simulated  trajectory,  and  (b)  tracking  error. 


Figure  11:  TL2  design  with  sensor  noise  s  and  the  disturbance 
nonlinearities,  (a)  Reference  and  simulated  trajectory,  and  ( 


(b) 


due  to  uncompensated  hysteresis  and  constitutive 
\  tracking  error. 


8.  Concluding  Remarks 

We  have  summarized  here  the  construction  of  macroscopic  constitutive  relations  and  models  to  characterize 
the  hysteretic  drive  dynamics  of  a  PZT-based  AFM  stage.  These  models  extend  previous  formulations  through 
the  use  of  general  density  representations  which  provide  high  accuracy  while  maintaining  implementational 
efficiency.  The  modeling  framework  also  facilitates  inversion  for  linear  control  design.  The  performance  of  initial 
robust  control  designs  for  high  accuracy  AFM  drive  regimes  is  illustrated  through  numerical  examples. 
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